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ABSTRACT 



Context. The puzzling detection of Can ions at fairly high latitude (i 30°) above the outer parts of the fi Pictoris circumstellar disk was 
, recently reported. Surprisingly, this detection does not extend to Na i atoms, in contradiction with our modelling of the emission lines in and 
nI ■ out of the mid-plane of the disk. 

Aims. We propose that the presence of these off-plane Can ions (and to a lesser extent Fei atoms), and the non-detection of off-plane Nai 
atoms, could be the consequence of the evaporation process of Falling Evaporating Bodies (FEBs), i.e., star-grazing planetesimals that 
evaporate in the immediate vicinity of the star. 

Methods. Our model is two-fold. Firstly, we show numerically and theoretically that in the star-grazing regime, the FEBs are subject to 
inclination oscillations up to 30 - 40°, and that most metallic species released during each FEB sublimation keep track of their initial orbital 
inclination while starting a free expansion away from the star, blown out by a strong radiation pressure. Secondly, the off-plane Can and Fei 
species must be stopped prior to their detection at rest with respect to the star, about 100 AU away. We revisit the role of energetic collisional 
processes, and we investigate the possible influence of magnetic interactions. 

Results. This dynamical process of inclination oscillations explains the presence of off-plane Can (and Fei). It also accounts for the absence 
of Nai because once released by the FEBs, these atoms are quickly photoionized and no longer undergo any significant radiation pressure. 
Our numerical simulations demonstrate that the deceleration of metallic ions can be achieved very efficiently if the ions encounter a dilute 
neutral gaseous medium. The required Hi column density is reduced to ~ 10 17 em"" 2 , one order of magnitude below present detection limits. 
• i-H . We also investigate the possibility that the ions are slowed down magnetically. While the sole action of a magnetic field of the order of 1 p.G is 
' not effective, the combined effect of magnetic and collisional deceleration processes lead to an additional lowering of the required Hi column 
S—j density by one order of magnitude. 
_ _ i 

Key words. Stars: circumstellar matter - Stars individual: Pic- Methods: analytical - Celestial mechanics - Methods: numerical - Planetary 
systems: protoplanetary disks - Molecular processes - Magnetic fields 
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1. Introduction 

The dusty and gaseous disk surrounding the young main- 
sequence star /} Pictoris has be en the subject of inte nse in- 
vestigation since its discovery (iSmith & Terrild Il984l) . The 
main motivation for these studies is that this disk constitutes 
the most convincing example of a probable young extrasolar 
planetary system, possibly analogous to the early solar sys- 
tem, p Pic is a young star, but it is not a pre-main sequence 
star. Its age has been subject to controversy in past years, 
but successive determinations base d on the kinematics of the 
socalled /3 Pictoris moving group dBarrado v Navascues et al 



1999tlZuckerman et al.l200ll:IOrtega et al.l2004 lead to a most 
recent estimate of 11.2Myr dOrtega et al.ll2.004l) . This shows 
that its disk should be called a young planetary rather than 
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a protoplanetary disk, meaning that planet formation already 
should have had enough time to occur. Indeed, although no 
direct planet detection has been made so far, several indirect 
observational facts suggest that planets are present in the disk. 
This mainly concerns asymmetries found in t he disk images 
from both scattered stel lar light by the dust dKalas & Jewitt 
1995b iHeap et al.1 12000) and thermal emission by the dust 



(Weinberger et al 



2003), that have been modelled as resulting 



from the gravitationa l perturbations by one Jupit er-sized plan et 
dMouillet et al.lll997tlHe"ap et alfeOOotlAugereau et alfeOOll) . 

The gaseous counterpart of the dust disk was det ected 
in absorption in the stellar spectrum dHobbs et all 1 1985b . and 
has been regularly observed since that time. Observations 
of many meta llic species such as Na i, Ca n, Fe n. . . have 
been reported (IVidal-Ma djar et al. 1994 ). e xtending more re 



cently to more fragile species like CO (Joll y et al 
Lecavelier des Etangs et al ] |200ll) . 
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The gas was first detected at rest with respect to the star, but 
Doppler-shifted, highly time-variable components are regularly 



observed in the spectral lines of many elements dFerlet et. al 



19871: lLagrange et all 1 19961: iPetterson & Tobinl 119991) . These 



transient spectral events have been successfully modelled as 
resulting from the sublimation of numerous star-grazing plan- 
etesimals (sev eral hundred per year) in the immediate vicinity 
of the star ( see iBeust et all 1 99811 1 996i and refs. therein). This 
scenario has been termed the Falling Evaporating Bodies sce- 
nario (FEBs). 

From a dynamical point of view, the origin of these nu- 
merous star-grazers seem to be related to mean-motion res- 
onances (mainly 4:1 and 3:1) with a Jovian planet orbit- 
ing t he star on a moderately eccen tric orbit (e' =s 0.07- 
0.1) dBeust & Morbidellill 19961 12000|) . In this context though, 
the suspected re sonance reservoirs are ex pected to clear out 
very quickly. In IThebault & Beusti (1200 II) . it was shown that 
collisions among the population of the planetesimals con- 
stituting the disk could help replenish the resonances from 
adjacent regions and subsequently sustain the FEB activity. 
However, there are still unsolved questions concerning this sce- 
nario. The main one concerns the amount of material avail- 



able ( Thebault et al. 20031) . There is a large discrepancy be- 
tween the number of planetesimals deduced from the FEB 
model and that deduced from an extrapolation of the dust ob- 
served pop ulation up to ki l ometr e-sized bodies. We neverthe- 
le ss note in IThebault et al J (12003) that the mass determination 
of IThebault & Beusti d200lh from the FEB scenario is very im- 
precise as it is indirect. Indeed if the collisions are supple- 
mented by some more violent transport processes, then the 
planetesimals population required to sustain the FEB activity 
could be much lower. 

An important outcome of this model is that it implies an 
important reservoir of planetesimals in the disk and the pres- 
ence of at least one giant planet at ~ 10 AU from the star. This 
is another argument in favour of the presence of planets in the 
B Pic disk. Moreover, the presence of numerous planetesimals 
is also a requirement of the dust models. Due to an intense radi- 
ation pressure, many dust particles should be quickly removed 
from the system. The particles observed consist of second gen- 
eration material continuously replenished from inside the disk 



by pl anetesimals, either by slow evaporatio n dLecavelier et al 



api 

19961) or by collisions ( Artvmowiczl 1997 ). This justifies the 



name second generation or debris disks given to the B Pic 
disk and other similar disks, s uch as the HD 141569 disk (see 
Augereau & Papaloizoul l2004, and refs. therein). 



Radiation pressure affects not only dust particles, but also 
the metallic species seen in absorption with respect to the star. 
Many of them undergo a radiation fo rce from the star tha t 



largely overcomes the stellar gravity (lLagrange et al.1 1 1 998 ). 



This is for instance the case of Can for which the radiation 
pressure is 35 times larger than the stellar gravity. This seems 
in contradiction with th e detection of ci r cumst ellar gas at rest 
with respect to the star. lLagrange et al. I d 19981) suggested that 
this stable gas was produced from inside by the FEBs them- 
selves, that it is then blown away by the radiation pressure, and 
afterwards slowed down by a dense enough H i ring where it ac- 
cumulates. The exact shape of this ring is of little importance, 



the main parameter being the integrated H i column density, of 
the order of 10 18 crrT 2 . Detailed modelling shows that all stable 
circumstellar lines can be reproduced this way. 



In a recent paper. lBrandeker et al.l (12004. hereafter B04I) re- 
port the detection with VLT/UVES of emission lines of metals 
(Fei, Nai, Can) in the B Pic disk, i.e., away from the direc- 
tion of the star. They report the detection of Na i and Fe i up to 
more than 300 AU fr om the star. Nai wa s reso lved earlier by 
Olofsson et al. (2001), jut the detection of B04 extends further 
out. H ? was also claimed to be detected in emission by lThi etakl 
(1200 ll) . implying huge quant ities (~ 50 M ffi ) in the /?Pic system , 
but this was questioned bv lLecavelier des Etangs et al.l (1200 lb . 
who reported from FUSE/LYMAN observations an upper limit 
Af(H_2) ^ 10 18 crrT 2 for the FL. column density towards B Pic. 



A particularly puzzling outcome of the B04 observations is 
the detection of Ca n emission at fairly high latitude above the 
mid-plane of the disk. Ca n is detected at 77 AU height above 
and below the mid-plane at 116 AU from the star. This corre- 
sponds to an inclination of 33° above the mid-plane. At this 
distance, Ca n is detected in both branches of the disk and on 
both sides of the mid-plane, and the emission at 33° inclination 
largely overcomes that in the mid-plane. Surprisingly, the Fe i 
and Nai emission do not exhibit such a structuring. It is con- 
versely concentrated in the mid-plane of the disk. However, the 
Fe i emission is broader than the Na i. At the height above the 
mid-plane corresponding to the peak emission in the Ca n lines, 
the Fe i is still detected. The gas shares the same radial velocity 
as the star within 1 or 2kms~ 1 at most. 

There is no straightforward explanation for the presence 
of species like Ca n at such latitudes above the disk, nor for 
the absence of other species. The purpose of this paper is to 
propose that this gas could constitute material released by the 
FEBs in the vicinity of the star, first blown away by radiation 
pressure, and then stopped far away from the star by friction 
with some gaseous medium, and/or by magnetic interaction. 
The Ca n and possibly the Fe i reach a significant inclination 
because the parent bodies (the FEBs) initially orbiting within 
the plane of the disk undergo inclination oscillations up so sev- 
eral tens of degrees when they reach the star-grazer state. Once 
released by the FEBs, the ions keep track of that inclination. 
In Sect. 2, we model the formation of the emission lines in and 
out of the mid-plane of the disk. We show that all the emissions 
are compatible with solar relative abundances between the el- 
ements under consideration, except that sodium should neces- 
sarily not be present (or strongly depleted) in the high latitude 
gas. In Sect. 3, we expose the dynamical model for high lati- 
tude gas generation, and we detail the theoretical background 
for inclination oscillations in the FEB state. We show that due 
to a negligible radiation pressure on Na n, sodium should not 
be present in this gas, in agreement with the observations. 

In Sect. 4, we investigate how the Ca n ions could be slowed 
down at the stellar distance they are observed. We show that 
this deceleration can be achieved by collision with a dilute neu- 
tral medium, and we discuss the role of elastic and inelastic 
collisional processes. We also discuss the effect of a non-radial 
magnetic field. While the magnetic field in itself is inefficient 
to decelerate the ions, we show that its presence increases the 
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Table 1. Measured circumstellar ratios, and simulated emission 
intensities for the transitions under consideration. The circum- 
stellar factors are taken from calibrated ESO/HARPS spectra of 
/3Pic (Galland F., private communication). The emission inten- 
sities (given in arbitrary units) out of the mid-plane are derived 
from the Cloudy run and those in the midplane are estimated 
by multiplying the former ones by the circumstellar factors. 



latitude), we may estimate the incoming kinetic energy flux F 
at 100 AU as 



Transition 


Circumstellar 


Emission 


Emission 




factor 


intensity 


intensity 






above the mid-plane 


in the mid-plane 


NaiD, 


0.91 


12.5 


11.4 


Nai D 2 


0.90 


24.8 


22.3 


CanK 


0.033 


15.8 


0.521 


CanH 


0.0057 


15.6 


0.0889 


Fe i 3859 A 


0.86 


23.9 


20.55 



efficiency of the deceleration by a neutral medium by an order 
of magnitude or more. 



2. Modelling the emission lines 

Our first task is to model the emission lines of Ca n, Na i and 
Fe i observed in and out of the mid-plane of the disk, as ob- 
served by B04. We focus on the Ca n K and H lines, the Na i 
Di and D2 lines, and the Fe 1 A — 3859 A line. Th e modelling is 
done u sing the radiative transfer code Cloudy bv lFerland et al.1 
d 19981) . 

We take the synthetic Atlas9 (Kurucz 199ll) stellar model 
for ft Pic with T eS = 8100 K, log g = 4.207 and a total luminos- 
ity of 8.7 L . We put gas at 1 16 AU from the star and compute 
the line emission. The chemical composition of the gas is as- 
sumed to be solar. The gas is modelled as a 10 AU wide layer 
with solar composition and a given hydrogen density. 

Without any additional energy source other than the stellar 
radiation flux, the emission in all lines as computed by Cloudy 
appears negligible, because the gas remains very cold (a few 
Kelvins). Thus, in order to generate detectable lines, the gas 
must be heated by some energy source. There can be some tur- 
bulence, but in the framework of our model, the strong decel- 
eration of the weak flux of incoming metallic ions might con- 
stitute a sufficient heating source. In the following we thus as- 
sume the emission lines to be excited thermally. The incoming 
ions, once blown away by the radiation pressure, reach a veloc- 
ity of ~ 10 3 kms _1 at 1 00AU. The variou s surveys of the FEB 
activity towards /3 Pic feeust et al. 1996b led to estimate that 
one roughly 1 kilometer-sized body is destroyed in front of the 
line of sight every day, with temporal fluctuations around one 
order of magnitude. Assuming that not all FEBs cross the line 
of sight, we estimate the total number of FEBs evaporated as 
N =~ 10 per day. Taking kilometer-sized bodies with an av- 
erage density of 2gcirT 3 , assuming that half of their mass is 
made of metallic ions that are pushed away by radiation pres- 
sure, and assuming that the ions spread over an open cone of 
a =~ 30° half-opening angle (in order to disperse ions at that 



F = -- 



1 Nmv 2 



0.2 erg cm 2 s 1 



(1) 



2 And 1 sin a 

where v ^ 1000 km s~' is the velocity of the incoming ions and 
m is the average mass of the FEBs. This kinetic energy heats 
the local H 1 gas over a distance corresponding to the stopping 
path of the metallic ions. This depends on the density of the 
local gas, but we see in Sect. 4 that an H 1 column density of 
~ 10 17 cirT 2 is expected to be sufficient to stop the ions over a 
distance of a few AU. 

We thus decided to perform runs of Cloudy with an extra 
heating source, with a volume-heating rate (parameter hextra) 
corresponding to the incoming kinetic energy flux deposited 
over the stopping distance. We performed several runs for hy- 
drogen densities ranging between 10 6 and 10 I0 crrr 3 . The re- 
sult listed in Table Q] are for 10 7 cm~ 3 . For other values, the 
absolute values of the emissions changes, but their relative be- 
haviour remains comparable. 

In almost all runs the same line behaviour is reported 
(Table [TJ: The Can K and H emissions are comparable, the 
Ca n K emission being slightly stronger; the Na 1 D2 emission 
is typically twice as strong as the Ca 11 K one, while the Fe 11 
emission is as strong as the Nai one. There is only little varia- 
tion between the "hot" and "cool" cases. 

This simulation is intended to hold for the gas out of 
the disk mid-plane. The spectrum of f3 Pic as seen from 
Earth is known to pr esent stable circumstellar components 



(ILagrange et al.l 1 19981) due to the gaseous counterpart of the 
disk. These componenents appear as sharp (~ 0.1 A wide) ad- 
ditional absorptions at the bottom of the rotationally broadened 
photospheric stellar lines. 

There is a major difference between ions located in and out 
of the mid-plane of the disk : the former ones see a stellar 
spectrum with these circumstellar components, while the lat- 
ter ones see a stellar spectrum without these components (be- 
cause they view the star across the disk, as Earth observers). 
The Atlas9 stellar model does not take into account these cir- 
cumstellar absorptions. Hence we must add them to the model. 
Unfortunately, adding additional spectral absorptions to the 
stellar models provided is not a standard procedure for Cloudy. 
It is thus not possible to accurately model the emission in the 
mid-plane of the disk. A first order approximation in order to 
take these absorptions into account is to apply reduction fac- 
tors of the stellar flux to ions in the mid-plane of the disk (and 
only to them). These factors are defined as the ratio of the stel- 
lar flux at the bottom of the circumstellar additional absorp- 
tions (if present) to the flux at the bottom of the corresponding 
photospheric lines (or equivalently the top of the circumstellar 
lines). They are simply measured from observed /? Pic spectra. 
We used ESO/HARPS spectra communicated by F. Galland. 
The factors are listed in Table Q] We note that the circumstellar 
lines are particularly deep for the Ca 11 lines. 

In order to derive a rough estimate of the line emissions, we 
may assume that the emission is proportional to the incoming 
flux. This only applies if the lines are optically thin, but the ratio 
of 2 between the Na 1 D2 and D; emission shows that this is the 
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case here. We apply the circumstellar factors of Table Q] to the 
emission intensities derived from Cloudy with no circumstel- 
lar absorptions. The ratio between Nai and Fei emissions re- 
mains unchanged (thanks to similar circumstellar factors), but 
the Ca ii emissions (in both lines) appears now far weaker (a 
few hundredth in relative intensity) than the Na i and Fe i ones. 

The observational constrains to fulfill (B04) are the follow- 
ing: In the mid-plane, the Can emission is small compared to 
the Nai and Fe i ones; the Nai D2 to Nai Di ratio is close to 2, 
showing that the emission is optically thin; the emission in the 
Fei A = 3859 A line is comparable to that in the Nai D2 line. 
Out of the mid-plane, the Ca 11 emission dominates, but the Fe 1 
line may still be detected, because the wing of the line is much 
larger than for the Na 1 lines; the Ca 11 K to Ca n H ratio is close 
to one. 

Our simple model succeeds in reproducing the emissions 
in the mid-plane. Out of the mid-plane though, the Nai emis- 
sion is always stronger than the Ca 11 one and comparable to the 
Fe 1 one. It is thus impossible to simultaneously have a strong 
Ca 11 emission, a Fe 1 emission a few times weaker, and an un- 
detectable Na 1 emission. The only possibility is to exclude the 
hypothesis of solar composition. 

We come therefore to the following conclusions : 1 ) for the 
emission in the mid-plane, the observations can be reproduced 
assuming solar relative abundances between iron, sodium and 
calcium; 2) out of the mid-plane, the Fe 1 and Ca 11 line intensi- 
ties can be consistently simulated assuming solar relative abun- 
dances; 3) the non-detection of Nai emission out of the mid- 
plane cannot be explained in these conditions, unless sodium is 
strongly depleted with respect to solar abundance. The model 
presented in next Section provides a plausible explanation for 
such sodium depletion. 



3. Falling Evaporating Bodies and high latitude 
ions 

3.1. General features 

We p ropose an origin for the high latitude ions observed by 
B04. Two facts need to be explained: i) Why are there large 
amounts of gas at 30° above the mid-plane of the disk ? ii) 
Why does sodium seem not to be present in this gas ? We pro- 
pose that this high latitude gas could be produced by the Falling 
Evaporating Bodies (FEBs). 

The FEBs are star-grazing bodies that fully evaporate in the 
vicinity of the star. Dynamically speaking, they are planetesi- 
mals that have been extracted from the disk orbiting the star and 
driven to high eccentricity orbits. They enter the FEB regime 
when their periastron reaches a threshold value (~ 0.4 AU) that 
allows the refractory material to evaporate. The details of the 
evaporation process of the bodies as their peria stron decreases 
down to a few stellar radii are described in iKarmann et al 
(2003). The bodies start to evaporate at each periastron pas- 



Star-grazers may be produced from a disk of planetesimals 
by planetary perturbati ons. The most efficient mechanism is 
the Kozai resonance ( Kozailfl962 ). which concerns bodies that 
have initial high inclination with respect to the orbital plane of 
the planetary system. Under secular perturbations, the initially 
highly inclined body is periodically driven to low inclination, 
but very eccentric, star-grazing orbits. This mechanism is re- 
sponsible for most of the sun-grazing bodies reported in the 
Solar System, such as the Kreutz group (Bailev et al. I ll992l) . 

However, the Kozai resonance is due to the secular, circu- 
lar part of the interaction Hamiltonian with the planet(s). It is 
therefore invariant with respect to any rotation in the planet's 
orbital plane. Bodies driven by the Kozai mechanism are thus 
expected to reach the FEB state with random orbit orientations. 
This does not match the statistics of the Doppler velocities of 
the variable spectral events observed, which shows a strong 
bias towards redshifts. Most of the suspected FEBs share some 
kind of common preferred periastron orien tation range whic h 
is not c ompatibl e with the Koz ai resonance dBeust etalJll996l) . 

iBeust & Morbidellil dl996l) proposed that the FEBs could 
be generated by another mechanism involving mean-motion 
resonances with at least one major perturbing planet. The secu- 
lar motion of bodies trapped in a given mean-motion resonance 
with a planet is usually characterised by coupled oscillations 
of the semi-major axis and the eccentricity around a median 
value, but if the planet's orbit is slightly eccentric, these os- 
cillations are superimposed on a long-term drift of the eccen- 
tricity that can in some cases bring it to star-grazing values. 



Yoshikawa ( 1989b showed that these chan ges are particularly 
important for resonances 4:1, 3:1 and 5:2. lBeust & M orbidelli 



( 1996) showed that the 4:1 reson ance is a potential source of 
FEBs via this mechanism, and in Thebaul t & Beust ( 2001 ). it 
was shown that the 3:1 may also contribute to the FEB phe- 
nomenon. The planet's eccentricity e' does not need to be 
very high; e' ^ 0.05 is enough, e' = 0.07 or 0.1 being typ- 
ical convenient values. Such eccentricity values are regularly 
reached by Jupiter due to its secular evolution. This mech- 
anism is close to the one that gave birth to the Kirkwood 
gaps in the asteroid belt, even if in the latter case, the over- 
lapping of mean-motion resonances wit h secular resonances 



considerably enhances the mechanism ([M orbidelli & Moons 

11995] 



sage, and their evaporation rate increases as the periastron dis- 
tance gets smaller. For the sizes considered (~ 10 km), the 
FEBs are fully evaporated when they reach a periastron value 
q ^ 0.15 + 0.05 AU. 



1993b iMoons & Morbidellil 1 19951; iMorbidelli & Moonsl 
Farinella et alJll994l) . We cannot exclude this enhancement as 
being effective in the f3 Pic system (this would imply the pres- 
ence of more than one planet), but there is no way to constrain 
it. We are perhaps witnessing in the /3 Pic system a process 
similar to what occurred in the Solar System, as it was of com- 
parable age to j6 Pics current age. 

Due to this mean-motion resonance mechanism, a given 
body, initially orbiting the star in 3:1 or 4:1 mean-motion res- 
onance with a Jupiter-sized planet, may reach the FEB state 
within ~ 10 4 orbital periods of the planet. Contrary to the 
Kozai case, the periastron longitudes of the FEBs at high 
eccentricity are now constrained by that of the perturbing 
planet, and share some common orientation in closer agree- 
ment with the observations. This scenari o was numerically 
tested over a large number o f particles dBeust & Morbi delli 
2000; iThebault & Beuslll2001l) . using the popular symplectic 



Beust and Valiron: High latitude gas in the f} Pictoris system 



5 



0.5 - 



10 2 10" 

time (years) 



3 10" 



40 



30 



20 



10 




10 2 10 

time (years) 



3 10" 



Fig. 1. Temporal evolution of the eccentricity (top) and of the 
inclination (bottom) of a typical particle trapped in 4:1 reso- 
nance with a planet orbiting p Pic at 10 AU with eccentricity 
e' = 0.07. The planet's mass is 1/1000 of that of the star. The 
initial eccentricity of the particle is 0.05 and its initial inclina- 
tion is 2°. 



integration package SWIFT_MVS (Wis dom & Holmanlll991 



Levison & Duncan! 119941) . It was shown that the suspected 



mechanism was able to fairly well match the statistics of the 
observed FEB velocities, provided the orbit of the perturbing 
planet adopts a given longitude of periastron with respect to the 
line of sight. If the disk of planetesimals holds a large enough 
population of bodies, collisi ons may help refill the r esonance 
and sustain the FEB activity dThebault & BeusfcoOl ). 

The orbits of the FEB progenitors in the mean-motion res- 
onances are supposed to be roug hly coplanar with the plan e of 
the disk. In the simulatio ns of iBeust & Morbidellil yOOO) and 
Thebault & Beustl ( 2001 ), the initial inclinations of the particles 



with respect to the orbit of the perturbing planet were initially 
chosen as less than 5° and 3° respectively, in order to mimic 
the typical distribution within a cold planetesimal disk. During 
their evolution within the resonance, as long as their eccentric- 
ity grows, the inclination of the particles remains small, but as 
they reach the FEB state close to e =s 1, their inclination is sub- 
ject to oscillations of larger amplitude, up to several tens of de- 
grees. This is illustrated in Fig.Q] which shows the secular evo- 
lution of the eccentricity and of the inclination of a typical par- 



ticle trapped in 4: 1 resonance with a planet orbiting /3 Pic. The 
particle starts at eccentricity e = 0.05 and evolves towards the 
FEB state at e =t 1 (the peak eccentricity is about 0.998), and 
then starts a decrease of its eccentricity. Of course the decrease 
phase is purely fictitious, as a real FEB would be destroyed 
by the successive periastron passages at peak eccentricity. The 
inclination, initially set at 2°, remains small for 10 5 yr and 
then starts oscillations that brings it far above the initial value, 
up to 40°. This oscillation regime stops after ~ 2.5 x 10 5 yr; 
it corresponds to the high eccentricity phase, when e ^ 0.85. 
Hence the FEBs may remain within the disk during most of 
their secular evolution, but finish in the FEB state with incli- 
nations that might bring them significantly out of the disk. The 
FEB themselves do not have excursions far out of the plane 
of the disk, because, as we will see below, when their inclina- 
tion is high, their argument of periastron a> is close to 0° or 
180°. This causes the major axis of their orbit to lie roughly in 
the plane of the disk at the time the inclination is high. As the 
orbit is very eccentric, the vertical excursion of the particle is 
limited. This is why this e ffec t, present in the simulati ons of 
Beust & Morbidellil J2000l) and lThebault & Beustl d200ll) . does 



not have much influence on the visibility of the FEBs (they 
need to cross the line of sight to be detected in absorption). 

What is true for the parent body is not necessary true for 
its byproducts. The metallic ions released by the FEBs such as 
Ca ii start to expand radially around the nucleus and stay for 
a while in a surrounding cloud that enables the FEB to be de- 
tected in absorption when it crosses the line of sight; but they 
are quickly expelled from there by the intense radiation pres- 
sure they suffer and then start a free expansion out of the sys- 
tem. This process is extensively descri bed in dedicated simula- 



Beust et al.l J 199dll 199611 1998b . Dynamically speaking 



tions in 

the ratio of the radiation pressure to stellar gravity is a constant 
(usually noted /J) for a given ion or dust grain. This is equiv- 
alent to the view that with radiation pressure, the ion feels the 
gravity of star as if its mass was multiplied by 1 - For Ca n, 
(3 > 1 (J3 — 35), so that the ions feel a negative mass and are 
strongly repelled by the star. They nevertheless follow a purely 
Keplerian orbit, namely a hyperbolic repulsive one, similar to 
the relative motion of two charged particles with charges of the 
same sign. This orbit is very different to that of the parent body. 
Both orbits share however the same orbital plane. Indeed, the 
ejection velocity of the mater ial escaped from the FEB (typi- 
cally ~ lkms -1 ; see refs. in IBeust et al. 1996) is very small 
compared to the orbital velocity of the FEB itself at a few stel- 
lar radii from the star (typically several hundreds of kms 1 ). 
The Ca n ions may thus be considered with reasonable accu- 
racy to have the same orbital velocity at ejection time as their 
parent body. Both Keplerian orbits are very different because 
of radiation pressure, but they share roughly the same orbital 
plane. In particular, the Can ions keep memory of the orbital 
inclination of their parent body at ejection time, even if this in- 
clination is large thanks to inclination oscillations in the FEB 
state. But contrary to the parent body, the orbit of Ca n ions are 
not confined close to the plane. As explained above, the argu- 
ment of periastron a> of the parent bodies is close to or 1 80° in 
the high inclination state; this forces the FEB to remain close to 
the mid-plane of the disk. Due to radiation pressure, the shape 



6 



Beust and Valiron: High latitude gas in the /3 Pictoris system 



of the Ca n ions is very different from that of their parent body, 
and their argument of periastron is not constrained in the same 
way. If they have a high initial inclination, they may evolve 
far off the plane of the disk as they escape from the system. 
If some dense medium is present at a given distance to brake 
them, they may be detected as an extended emission signifi- 
cantly above the plane such as in the B04 observation. Their 
detection of Can at 33° off the mid-plane could then well cor- 
respond to ions that have been produced close to the star by 
FEBs with similar inclination, and that have freely escaped up 
to 100 AU before being stopped there. 

Then, why should this process only con cern C a n and not 
the other species detected in emission by |B04| ? Iron and 
sodium are byproducts of dust sublimation in FEBs like cal- 
cium. Contrary to Ca n, Na i and Fe i are quickly photoionized 
by the star in the FEB environment. Fen (like Fei) still un- 
derg oes a radiation pres sure that overcomes the stellar grav- 



ity ( Lagrange et alJI 19981) . so that iron is expected to behave 



like calcium. However, unless the electronic density is high, as 
is the case in the vicinity of the mid-plane disk, iron remains 
predominantly i n the Fe n state whose spectral lines were not 
searched for bvlBOl. Conversely, Nan does not feel any no- 
ticeable radiation pressure, so that once produced, the Na n ions 
keep following the original orbit of the FEB. They may after- 
wards diffuse slowly in the mid-plane of the disk, but they are 
not subject to a quick off-plane ejection like the other species. 
In this context, we thus expect the gas expelled off-plane by this 
process to contain calcium and iron with solar relative abun- 
dances, but no sodium. This matches our analysis of the emis- 
sion lines. 

In the following we detail the theoretical background of the 
origin of the inclination oscillations of the FEBs in high ec- 
centr icity regime, and we show exa mples from the simulation s 
from lBeust & Morbidellil (l2000h and lThebault & Beustl J200lh . 

3.2. Three-dimensional motion in mean-motion 
resonance 

The theoretical background for t he FEB dynamics in mean - 
motion resonance is described in Beust & Morbidellil dl996h . 
but restricted to the planar case. Here we wish to extend it 
to three-dimensional motion. We describe the restricted three- 
body problem, i.e., a problem where a mass-less test particle 
orbiting a star is perturbed by a planet orbiting the star on an 
unperturbed Keplerian orbit. 

The full analytical analysis is presented in Appendix A. 
We assume that the particle is locked in a (p + q) : p mean- 
motion resonance with the planet. The resonant motion is 
usually described by the "c ritical angle of the resonance" cr 
dMoons & Morbidellill 19951) . with 



P + <1 y P , 

cr — A A ■ 

q q 



(2) 



where A is the mean longitude of the particle along its orbit; A' 
is the same for the planet; m is the longitude of periastron. 

Non-resonant orbits are characterised by a more or less reg- 
ular circulation of cr, while resonant orbits exhibit libration 



of cr around a stable position. If the planet's orbit is circu- 
lar, then the following quantity is a secular constant of motion 
dMorbidelli & Moonsli 1 9931; iMoons & Morbidellill 1 995b : 



- Vl -e 2 



COS I 



(3) 



where e is the eccentricity, i is the inclination, and p. is the mass 
parameter, i.e. the ratio of the planet's mass to the total mass. 
The inclination oscillations are controlled by this parameter. If 
the planet's orbit is not circular, then strictly speaking N is no 
longer constant, but as the planet's eccentricity <?' is moderate, 
its variation is slow. 

Considering the circular problem is equivalent to expand- 
ing the Hamiltonian 'H in powers of e' , and retaining only the 
leading term. The leading (circular) term is responsible for the 
cr-libration, but also for the inclination oscillations at high ec- 
centricity. The higher order terms in <?' cause a slow drift of the 
N parameter. This drift can drive the particle to high eccentric- 
ity, even in the planar problem; this is the way the FEBs are 
generated. 

Finally, the dynamics of the particle is characterised by 
three time-scales : a first, small one related to the cr-libration; 
a second, larger one characterising the inclination oscillations; 
a third, long one describing the secular eccentricity changes 
from ~ to ~ 1 . The second time-scale is larger than the first 
one, but significantly smaller than the third one. Hence during 
one inclination oscillation, the value of N may be considered 
as ~ constant, which is equivalent to considering the circular 
problem. 

In the circular problem, the secular Hamiltonian H depends 
only on the inclination ; and on the argument of periastron to, 
once the value of N is fixed. It is possible to explore the dy- 
namics ju st drawing level curves of fi in (w, plane, exactly 
as done in lBeust & Morbidellil (1 19961) . This is done in Fig.[2]for 
the 4: 1 resonance, for 4 different values of N. Instead of giving 
the value of N, we give a value for the eccentricity and com- 
pute the value of N that gives this eccentricity value for i - 0. 
The value of p was fixed to 0.001 as typical for a Jupiter-sized 
planet. The value of N is indicated above each plot, and a cor- 
responding eccentricity scale is given to the right of the plots. 
The four plots corresponds to eccentricity values at i — of 0.5, 
0.8, 0.9 and 0.99 (these values appear at the lower right corner 
of the plots). In all these plots, the eccentricity variation over 
most of the plot is very moderate. Hence each of these plots 
should be regarded as a picture of the dynamics in a given ec- 
centricity r egime. The plots of Fi g. [2] are equivalent to those of 
Fig. 5 from lMorbidelli & Moonsl (1 1993b for the 2:1 resonance. 

Consider now a given particle trapped in the 4:1 reso- 
nance which starts its eccentricity growth. In a low eccentricity 
regime, like the one for e - 0.5 (N = 1.973), the inclination 
remains low while a> circulates. Following a level curve of 
if the inclination is initially low (a few degrees), it undergoes 
small variations that keep it in the same range. At a higher ec- 
centricity regime, the phase portrait changes. We note in Fig. [2] 
that two islands of libration for co appear around to = ±90°. 
However, these islands of libration do not concern the particles 
we are considering. Our particles start within the plane of the 
disk with an inclination that does not exceed a few degrees. 



N=1.973 



Beust and Valiron: High latitude gas in the f} Pictoris system 

4:1 N=2.14 



4:1 



20 ; 



o 

CD 
Jh 
Oil 
CD 
13 



10 




SO. 3881 



0.4761 



100 200 300 

co (degrees) 
N = 2.244 4:1 



0.5 




100 200 300 

co (degrees) 

N = 2.429 4:1 



60 i 



m 
o 

0) 
Of) 

(1) 



40 



20 




SO. 4899 



0.7349 



0.8223 



80 



60 



0.8641 ^ 40 



0.8859 



0.8967 



20 



100 200 
co (degrees) 



300 



0.9 



t 
e 



100 200 

co (degrees) 



300 



0.3587 



0.6217 



^0.7211 



-0.7696 



-0.793 




0.99 



t 
e 



Fig. 2. Level curves of the Hamiltonian fi in the (to, i) plane for the circular problem, for particles having negligible cr-libration 
amplitude, for different values of the constant parameter Af. The eccentricity scale (denoted "e" on the right of the plots) is related 
to the inclination scale to the left by the constant value of N. 



Hence the curves they follow are those located below the is- 
lands of libration. For our particles, to still circulates, but fol- 
lowing the level curves, the inclination i is subject to periodic 
jumps up to possibly several tens of degrees when to reaches 
or n. The higher the eccentricity regime, the higher the incli- 
nation jumps. This is the origin of the inclination oscillations 
reported in the numerical integration. 

This dynamics is a resonant version of the Kozai dy- 
namics. In the non-resonant circular restricted problem, the 
Kozai Hamiltonian describes the secular dynamics of the par- 
ticle. It is obtained by a double averaging of the original 



Hamiltonian over the orbital motions of the planet and of 
the particles (IKinoshita & Nakailll999l) . It is well known that 
this Hamiltonian has a secular constant of motion which is 
the z-component of the angular momentum (or equivalently 
Vl - e 1 cos i = est). It is also well known that at high incli- 
nation, this Hamiltonian shows two islands of libration in (/, to) 
space around to = +7r/2, and that particles moving in these is- 
lands evolve periodically from a high inclination and low ec- 
centricity regime to a low inclination and hig h ecce ntricity. 
This behaviour constitutes the Kozai resonance (Kozaii il962l) . 
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Beust & Morbidellil d2000b and lThebault & Beustl d200ll) . This 
is illustrated in Fig. [3] I n this figure, we consid er a typical sim- 
ulation described in lBeust & Morbidellil d2000l) . with e' = 0.07 
and fi = 0.002. As described in that paper, the simulation is 
made over 10 4 particles initially chosen orbiting the star in 4: 1 
mean motion resonances with the perturbing planets. The ini- 
tial eccentriciti es are randomly chosen < 0.1 and the inclina- 
tions < 5°. In iBeust & Morbidellil d2000h . it was shown that 
many of the particles evolve to the FEB regime, and we per- 
form statistics over the expected Doppler velocities when the 
FEBs cross the line of sight. This statistic appears in agreement 
with the observational one. In Fig. [3] we display information 
about the inclinations of the particles when they are in the FEB 
regime. More specifically, we count all periastron passages for 
the particles with periastron values less than 0.4 AU (the FEB 
regime), and that are not yet destroyed by evaporation. Figure[3] 
shows that the highest inclinations correspond to the smallest 
periastron values. This is in agreement with the theory outlined 
above (Fig. |2} which shows that the high inclination oscilla- 
tions are to be expected in the high eccentricity regime only, 
i.e. in the FEB regime. Fig. [3] also shows that whenever they 
reach high inclinations, the FEBs assume an argument of peri- 
astron close to +7r/2. This is in agreement with the theory, and 
demonstrates the reality of the proposed mechanism. Moreover, 
the concentration of points close to a> ± n/2 shows that the 
FEBs spend more time in high inclination regimes than in low 
regimes. This is confirmed directly by the temporal evolution 
of the inclination in Fig.Q] 



Fig. 3. A statistical test of the inclination oscillation regime for 
FEBs. These results conce rn a typical simulation described in 
Beust & Morbidellil d2000l) with e' = 0.07 and fi = 0.002. Each 
dot corresponds to a series of periastron passages of a body 
that has entered the FEB regime (periastron ^ 0.4), and that is 
not yet fully evaporated. Top plot : Inclination as a function 
of periastron; Bottom plot : inclination as a function of the 
argument of periastron a> 



The islands of libration in the plots of Fig. [2] describe a 
Kozai resonance, within a mean-motion resonance. Indeed, as 
a is fixed the condition N = est is exactly equivalent to the 
Kozai condition e 2 cos 2 = est. The FEBs trapped in the 
4: 1 resonance that evolve at very high eccentricity regimes are 
concerned by this, but they are not trapped into the Kozai reso- 
nance, as their argument of periastron a> still circulates, and as 
they periodically return to i =s 0. However the Kozai dynam- 
ics influences them and causes periodic inclination jumps up to 
several tens of degrees, even if the initial inclination is low (a 
few degrees). 

3.3. Tests over a large number of bodies 

To test the statistical effect of the eccentricity jumps reported, 
a numerical test over a large number of bodies is necessary. In 
fact it was not necessary to perform new simulations. We just 
take the (still available) results of the simulations described in 



4. Stopping the ions 

We showed in the previous section that, due to their inclina- 
tion oscillations, the FEBs constitute a potential source of Ca n 
and Fe i ions that may escape far off the disk plane. These ions 
- being blown away by a strong radiation pressure - should 
accelerate mostly freely and reach the outer parts of the disk 
with high velocities, up to ~ 1000 km s~'. However, such ve- 
locities would be in sharp contradiction to the observations by 
B04, where the Can ions are found at rest with respect to the 
star, about 116 AU away. This implies that some process is able 
to efficiently slow down the ions despite the intense radiation 
pressure they undergo (J3 = 35, see above). 



This issue was investigated recently by Fernandez et al 



( 200 6. hereafter F06I) . as many metallic species are observed at 
rest relative to the star despite a strong radiation pressure. They 
identified three possible braking processes : collisions among 
ions, collisions with charged ions, and collisions with a neu 



tral g as. Collisions among ions are very efficient dBeust et al 



Il989h . and the whole plasma tends to behave like a single fluid 
with a weighted average [3. However. lF06l show that unless car- 
bon is overabundant, the fluid is still accelerated by the star. 
Collisions with charged grains are only efficient if the grains 
are mostly carbonaceous. Moreover, at high latitude in the disk 
where the Ca n is observed, the density of the dust is low (0.3% 
of that of the midplane according to the profile given by|F06J). 
Coll isions with neutral gas are conversely a good candidate. 
IF06i showed that a minimum mass of neutral gas of ~ 0.03 M® 
is enough to stop the ions. However, due to the high incom- 
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ing velocity of the Can ions, the basic analytic model must 
be revised. We detail this below, and conclude that the actual 
braking is even more effective than in the basic formulation. 



4.1. Stopping with gas 



Neutral species like H i or He i or even H2 ar e not subject to 



any si gnificant radiation pressure from the star dLa grange et al 



1998); they may thus stay orbiting in a Keplerian way at some 
distance from the star. The incoming Ca 11 ions may then col- 
lide into this buffer gas and be slowed down to negligible ve- 
locity. This process was also invoked as a way to stop Can 
ions a few AU from the sta r in order to generate the stable cir- 
cumstellar Ca 11 absorption dLa grange et al. Il998h . In that work, 



we showed that a column density of ~ 10 cm is enough to 
slow down the Ca 11 ions. Here we investigate whether the same 
mechanism could also apply to decelerate faster off-plane ions. 



The rapid off-plane ions are not stopped at a few AU like 
those that stay within the plane because they do not encounter 
any noticeable gaseous medium at their orbital inclination. 
Why should they be stopped around 116 AU ? We must as- 
sume that at such a distance, the disk tends to flare. Hence 
some dilute material could be present at 30° or more in the 
outer disk while remaining absent in the inner disk. But why 
116 AU ? This distance corresponds approximately to the loca- 
tion of the power l aw break-up in th e surface brightness radial 



-up in the 
ai] l2000h . 



profile of the disk (Hea p et al.ll2000t) . Closer to this threshold, 
the surface brightness decreases as r _u , while further away it 
falls off much more s teeply as r~ . This was interpreted by 



Augereau et al.l d2001l) as a consequence of the distribution of 
planetesimals in the disk. The dust particles are produced by 
the planetesimals and then scattered into the outer disk by the 
stellar radiation pressure. The power law break-up at ~ 120 AU 
is consistent with a planetesimal disk presenting a rather sharp 
outer edge located at this distance dAugereau et al ]|200lb. The 
planetesimals disk appears thus to be truncated at the same dis- 
tance where the off-plane Ca n ions stop. These facts may be 
related. The flaring of the disk that we invoke at that distance 
for stopping the Ca 11 ions could be due to t he perturbations by 



succe ssive stellar flybys, as was invoked bv lLarwood & K alas 



(2001J) as an explanation for asymmetries and arc-like struc- 
tures in the outer parts of the circumstellar dust disk. But the 
same mechanism could also be invoked to account for the trun- 
cation of the planetesimal disk at the same distance. Further 
away than 120 AU, the planetesimals are perturbed by stellar 
flybys, and may not remain in a thin disk. Also the planetes- 
imals themselves may not have had the opportunity to form 
there. The stellar flybys may have scattered away (and proba- 
bly in the vertical direction) the initial material from which the 
planetesimals were expected to form. Thus, the /? Pic disk be- 
yond 120 AU could still be in a kind of primordial state where 
no refractory material condensation would have occurred, with 
a significant flaring due to stellar flybys. 



4.1 .1 . The analytical induced dipole model and its 
limitations 

We now review the basic mechanism of decelerating by a neu- 
tr al medium, noting that we depart from the situation described 
in lLagrange et alj dl998l) . in two points: i) the volume density 
of the incoming ions and of the colliding medium is proba- 
bly much less at 1 16 AU and 30° inclination than in the plane 
at a few AU; therefore no significant pressure effect is to be 
expected and the hydrodynamic description may be dropped; 
ii) the incoming velocity of the Ca n ions is likely to be much 
higher. In the limiting case of a free ballistic runaway driven by 
the radiative pressure of the unobscured star, the resulting ve- 
locity is ~ lOOOkms" 1 , about two orders of m agnitude higher 
than was considered in lLagrange et al.l d 1998D . 

A simpl e collisional decele rating mechanism was initially 
described in Beust et al. (1989) for moderate velocities. When 
a charged ion approaches a neutral atom, a dipole is induced 
on the neutral atom, from which an interaction results between 
the two particles that may be well described by the potential 
energy 



V(r) = 



1 



aq 



4ne 2r 4 



(4) 



where a is the polarizability of the atom, r is the interaction 
distance, and q is the ch arge of the ion, th e other symbols hav- 
ing their usual meaning dMcDanielll 1 9641) . The relative motion 
in this potential cannot be solved exactly, but there is a criti- 
cal impact parameter bo (that depends on the impact velocity 
v) separating two regimes : for b > bo there is a minimum 
approach distance between the particles preventing any close 
collision; for b < bo this is not so, and the particles undergo a 
physical collision (see appendix B). A fairly correct approxi- 
mation is then to neglect the effect of the encounter for b > bo, 
and to consider that the mean impulsion loss by the ion in the 
physical collision is fxv, where /i is the reduced mass. The in- 
teraction cross section is then nb\ and the resulting drag force 
/ is opposed to velocity : 



/ = -/jTib^nvv 



(5) 



where n is the number volume density of neutrals. The expres- 
sion of bo (see appendix B) is 

.2\l/4 



b 



1 A-aq 1 
4neo nv 2 



(6) 



so that the resulting force is proportional to the velocity: 
/ = -kv with 



k = nn 



4-jiaq 2 



\ 47T60 



(7) 



We confirmed to within 5% this simple analytical approxima- 
tion of the drag force using the more general numerical treat- 
ment described later on. 

However this induced dipole regime holds as long as the 
drift velocity v is not too large. When v grows, bo becomes 
smaller than the physical size of the particles. In this regime the 
interaction cannot be described any longer as attractive, and it 
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Table 2. Values of the critical impact parameter bo as a function 
of the relative velocity v, as computed from Eq. (O for the Ca n 
- H i interaction 



v (km s 1 ) 
bo (nm) 



1 

0.785 



10 

0.24 



100 
0.078 
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0.024 
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Fig. 4. The ab initio calculated interaction potential between 
Ca ii and H i (in their ground triplet state) as a function of the 
mutual distance (fat line), superimposed on the induced dipole 
potential (thin line) 



rather approaches a hard sphere regime at shorter range with a 
constant cross section that does not depend on v. According to 
Eq. (0, the drag force turns out now to be proportional to v 2 
instead of v. Table [2] lists computed values of bo for different 
values of v for the Ca n - H i interaction (the polarisability of 
Hi is 6.7xl0~ 31 m 3 ). We note that as soon as v ^ lOkms -1 , this 
induced dipole model becomes unrealistic because bo becomes 
comparable to or smaller than typical atomic radii. Hence Ca n 
ions that encounter H i atoms at ~ 1000 km s _1 undergo a colli - 
sional interaction that is similar to a hard sphere regime. As the 
velocity decreases along successive collisions, the interaction 
finally enters the induced dipole regime. A correct description 
of the decelerating process of the Ca n ions implies therefore to 
be able to describe the interaction at every velocity, in particu- 
lar in the intermediate velocity regime between the two above 
described extremes. 

4.1 .2. The "smooth sphere" model 

In order to have a more coherent description, we introduced a 
"smooth" sphere approximation based upon a continuous de- 
scription of the interaction potential U(r) between Can and 
Hi as a function of the relative distance r. The interaction 
originates from a quantum-mechanical interaction at the mi- 
croscopic level. During a collision, the incoming Ca n and H i 
particles form an intermediate molecular ion. This molecular 
complex possess two valence electrons originating from each 
incoming particle. In a quantum description of the interaction, 
these two electronic spins recouple to form either a singlet or 
triplet state with a 3 to 1 probability in favour of the triplet 
state. However in the case of energetic collisions these input 



states will interact with various excited states of the same spin 
multiplicity, leading to numerous inelastic processes that will 
be sketched in the the next subsection. 

Fortunately the ground state triplet state is not expected to 
be very reactive for moderate collisional energies and could 
provide a realistic basis for our "smoo th" sphere mode l. Using 
the ab-initio Gaussian 94 package ( Frisch etall 1 1 995b we in- 
vestigated the triplet input state using a restricted open-shell 
approach (in order to preserve the total spin). We performed 
ROMP2 calculations to take into account the electronic correla- 
tion and also to some extent the interactions with excited states. 
We added diffuse and polarization functions on both centers, 
paying special attention to the proper description of the polar- 
isation of the hydrogen atom. We also performed a population 
analysis at the Mulliken level to check that charge transfer ef- 
fects were small (contrary to the singlet state where strong in- 
verse charge transfer effects were found). The resulting triplet 
potential is plotted in Fig.|4]as a function of r, superimposed on 
the induced dipole interaction. Our ab-initio potential is consis- 
tent with the induced dipole approximation beyond 0.5 nm, but 
closer to 0.5 nm, it exhibits a repulsive wall providing a smooth 
transition towards the limiting hard sphere regime. 

Let us now take into account this "smooth" sphere poten- 
tial for a determination of the drag force. Given any interaction 
potential U(r), the drag force acting on the ion is opposed to 
the velocity and may be written as 



(8) 



where b is the impact parameter, and x(b, v) is the deflection 
angle due to the encounter corresponding to b and v. x(b, v) 
itself depends on the form of the potential and may be obtained 
from an integral along the relative motion during the encounter 
(see appendix B). 

4.1 .3. Beyond the "smooth sphere" model: the 
"inelastic" model 

The above expression for the drag force assumes implicitly that 
the encounters between the ion and the hydrogen atoms are 
elastic. 

On the contrary, energetic collisions are likely to trigger 
various inelastic processes. For impact velocities in the range 
100-1000 kms" 1 , the energy available in the center of mass is 
huge, from 50 eV to 5000 eV, and is able to induce a large vari- 
ety of excitations in both the valence and core electronic space. 
Beyond the mere electronic excitation of either particles, these 
energetic collisions can thus trigger a whole range of inelastic 
and reactive processes including inverse charge transfer, single 
or multiple electronic ionizations, etc. . . Collisions with He i or 
H2 will be even more energetic due to the larger reduced mass. 
A detailed theoretical description of all these processes and of 
their cross sections and branching ratios is beyond the range of 
the present study. Experimental investigations might provide a 
better starting point for further studies. 

All these inelastic processes will convert a part of the in- 
coming kinetic energy into internal energy of the particles and 
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also, if ionization occurs, into kinetic energy of the secondary 
electrons. Of course the internal excitations will mostly decay 
radiatively and will never be converted back to kinetic energy 
of the Can ions. In addition, the Can ions are likely to in- 
crease their charge, following either simple or multiple ioniza- 
tion, or inverse charge transfer with H iLj. Once ionized to Ca m 
or higher, the calcium ions are expected to recombine to Ca n 
after a while. But during the time they spend at higher ioniza- 
tion states, they should be decelerated even more effectively, 
because they no longer feel any significant radiation pressure 
from the star (as the species under consideration have no strong 
spectral lines in the visible-UV domain), while the drag force 
is expected to increase with ionization level. 

Moreover, the neutral H i gas, once shocked by the incom- 
ing high velocity ions, is expected to be partly ionized by this 
process. The gas will thus tend to behave like a plasma with 
a collective dynamical behaviour, resulting in an averaged [3 
ratio, as described by lF06l These collective effects should en- 
hance significantly the efficiency of the braking process. 

All these inelastic effects will increase the energy loss 
by the Can ions and thus the drag force predicted by the 
smooth sphere model. Their description is beyond the scope 
of the present paper. In the following we propose an order-of- 
magnitude calculation using a very crude model to describe the 
ionization processes involving collisions between Ca n ions and 
H i atoms. A simple way to treat possible ionization is to mon- 
itor the available kinetic energy 1 12jiv 2 before the collision in 
the inertial referential frame. First we select close collisions for 
which the interaction departs from the induced-dipole model. 
Second we consider the opening of successive ionization chan- 
nels when the energy is augmented. We model the ionization 
for H i and up to 7 electrons for Ca n, assuming an average en- 
ergy loss I p = 20 e V per electron. This loss is supposed to take 
into account both the extraction energy and the kinetic energy 
of the expelled electron. The ionization limit of 7 electrons for 
Can is rather arbitrary and includes the 2p shell and the outer 
3 s electron. The most energetic collisions might rather ionize 
an inner 1 s electron, but this would result in a comparable en- 
ergy loss because their binding energy is higher, about 150 eV. 

In practice, we modify the smooth sphere model with the 
following prescription. If the available kinetic energy exceeds 
k x I p with k < 8 and if the closest approach between the two 
particles is less than 0.5 nm, the incoming kinetic energy is ar- 
bitrarily reduced by kxlp. This causes the relative velocity after 
the encounter V to be less than the initial velocity v. We define 
the energy restitution coefficient e < 1 as v' — ev. This may be 
written as: 



kL 



1 a 1 2 2 1 2 
—mv — —e mv — —mv 

2 2 2 

If we assume for simplicity that the deflection angle x is un- 
changed with respect to the elastic case, then the expression of 
the force becomes now 



f — -Tjinyw 



. ( 



2e sin - 



+ 1 -e\bdb 



■(10) 



1 In view of our ab-initio investigations, this latter process is ex- 
pected already to be important at moderate energies for collisions in a 
singlet electronic state. 
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Fig. 5. The drag force on Can ions due to a Hi medium of 
unit density (1 irT 3 ), as a function of the drift velocity v. The 
thin line corresponds to the induced dipole model, the thick 
grey line to the smooth model (based upon the triplet potential 
shown in Fig. |4j and the thick black line to the inelastic model 
(based upon the combination of the triplet potential and a crude 
treatment for ionization effects, see text). 



Here the coefficient e is implicitly a function of b and v, as 
explained above. 

The result of the force computation in the various cases is 
shown in Fig. [5] as a function of the relative velocity v. As in 
any case the drag force is proportional to the density of the 
H i medium encountered, showing the force for a unit density 
medium is enough for comparison purposes. In the induced 
dipole approximation case, we find as expected / oc v, and we 
see that this approximation is valid up to v - 20 km s . At the 
higher velocity regime, we have / oc v 2 when the triplet poten- 
tial is taken into account, corresponding to our smooth sphere 
regime. When ionization is also taken into account, the non- 
elastic character of the interactions adds an extra force term 
to the elastic smooth sphere case. Inelastic effects turns out to 
be particularly noticeable for 100 kms"' $ v $ 1000 kms~' 
(at higher velocity ionization is present but the smooth sphere 
regime dominates). This velocity regime concerns Ca n ions en- 
countering a Hi medium at about 116AU. Consequently, this 
order-of-magnitude calculation suggests that a proper inclusion 
of ionization and other inelastic effects would significantly en- 
hance the effective drag force beyond the predictions of the 
smooth sphere model. 



( 9 ) 4.1 .4. Estimate of the required H i column density 

Now we have an estimate of the drag force, it is of interest to 
derive in which conditions the H i medium is able to stop the 
Can ions. Irrespective of the initial velocity, if the medium is 
dense enough, the ions will always be stopped. The question 
is to know the amount of H i neutrals required to do this, and 
more specifically the column density N s the ions need to cross 
before being stopped. Let us first consider a simplified case 
where radiation pressure is not taken into account. The ions 
arrive at initial velocity vq and encounter an H i medium with 
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volume density n. The equation of motion of an ion along its 
path will be 



dv 

m — = —nF (v) 
d? 



(11) 



where m is the mass of the ion, and F(v) is the force per unit 
density shown in Fig. [5] We change the dependent variable t to 
the integrated encountered column density N (dN = nvdf): 



dv 
cDv 



-F(v) 



(12) 



from which we immediately derive the integrated column den- 
sity N s required to slow the ions from vq to 0: 



N, 



r 

— m I 

Jo 



F(v) 



■dv 



(13) 



Let us now reintroduce the radiation pressure as a constant 
force P. The equation of motion is 



dv 

m — = P — nF(v) 
dt 

which is equivalent to 

dv P 

mv — = F(v) 

dN n 



(14) 



(15) 



The ions are now no longer exactly stopped, but rather slowed 
down to an equilibrium velocity v eq characterised by nF(v eq ) = 
P. In practice this terminal velocity is low, so that the ions may 
be considered as stopped. The radiation pressure P at 116 AU 
on Ca ii ions is ~ 1 .7 x 10~ 29 N (35 times stellar gravity). Let us 
take the uppe r limit of 10 18 cm" 2 for the hydroge n column den 



sity given bv lLecavelier des Etangs et"aL ( 2001 ). spread over a 
distance d. We derive a ratio 



- = ^(v eq ) - 2.5 10 
n 



-41 



x d(AU) N nT 



(16) 



If we consider as a maximum value for d a few tens of AU, a 
comparison with Fig.|5]shows that v eq is less than 1 kms" 1 and 
that it falls well within the induced dipole regime. This result 
still holds even if we assume a column density lower by one or 
even two orders of magnitudes. 

Strictly speaking, an infinite column density is required to 



reach v eq . We should write 



r 

— m I 



F(v) - F(v eq ) 



dv 



(17) 



As close to v eq , we have F(v) oc v, we derive that the integral 
diverges logarithmically towards v eq . But reaching a velocity 
that is of the same order of magnitude as v eq is enough for our 
purpose. Moreover, for v >> v eq we have F(v) » F(v eq ), so 
that F(v eq ) can be neglected in Eq. [17] Finally, Eq. ( fT3T > turns 
out to be a good estimate for N s even in the presence of radia- 
tion pressure. This is due to the fact that the terminal velocity 
v eq is very low. 

A^ as given from Eq. ( fT3l is plotted on Fig. [6] as a func- 
tion of the initial velocity, for the various interaction mod- 
els considered. As expected, in the induced dipole regime, 
we have N s oc vq, but at higher velocity, N s is reduced by 
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Fig. 6. Hi column density necessary to stop Can ions, as a 
function of the initial velocity Vo, as derived from Eq. Qj] ac- 
cording to several models (see text). The plotting conventions 
are the same as in Fig. [5] 

several orders of magnitude with respect to that crude esti- 
mate. The smooth sphere model causes N s to stay below a few 
10 17 cm" 2 (asymptotically N s oc In v ). With v = 1000 km s" 1 , 
we predict N s 10 17 cm" 2 . This is one order of magnitude be- 
low the upper limit to the H? co lumn density towards ft Pic 



(ILecavelier des Etangs et al.l 1200 lb . As suggested by our in 



elastic model, the inclusion of inelastic effects would further 
lower the require d colu mn density. Moreover, the collective ef- 
fects decribed by lF06l due to partial ionization of the neutral 
gas, are expected to enhance the braking process. N s could thus 
be even less than the value we derive. 

Hence we stress that the model we present here provides a 
plausible mechanism for stopping the Ca n ions at 100 AU from 
/3 Pic, in order to render them detectable in emission. 

4.2. Stopping with a magnetic field 

Another natural interaction ions that may be subjected to is 
the influence of a magnetic field B they would encounter at 
100 AU. In the presence of a magnetic field, the equation of 
motion of an ion is 



dv „ 
m — = P + qv X B 

dt 



(18) 



where P is the radial radiation pressure. Let us consider for 
simplicity that P is constant (the ions are supposed to be 
stopped over a short distance), as is B. The equation of motion 
has a well known analytical solution. Introducing a referential 
frame (Ox, Oy, Oz) where B \\ Oz and P lies in the xOz plane, 
and the angle <p between P and the xOy plane, the solution may 
be described in that referential as 



V v = 0) g COS ( 

Vy — OJg COS ( 



R cos(w^f) + L sin(w g f)] 
L(cos(w^f) — 1) - R sin(wyf)J 



(19) 



) (Lojgt + R} 



The ion is supposed to initially move radially at velocity vo. 
Here co g = qB/m is the gyromagnetic frequency; R = mvo/qB 
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is the gyromagnetic radius associated with the velocity vo; L 
is defined as P = mLto 2 , i.e., it is a characteristic length asso- 
ciated with the strength of the radiation pressure. The motion 
perpendicular to the field is a cycloid-like motion, i.e., a com- 
bination of a circular motion of radius p = V/? 2 + L 2 cos <p at 
frequency u> g and of a linear drift at velocity Lco g cos <p. The 
motion parallel to the field is uniformly accelerated by the ra- 
diation pressure. It is important to note that here, contrary to 
the neutral drift model, the velocity does not reach an asymp- 
totic value. If the field is not perpendicular to P (sin 0^0) 
the velocity even increases continuously. But even if sin <p — 
the motion in the xOy plane is still cycloid-like, and the mod- 
ulus of velocity undergoes a periodic modulation of amplitude 
2pio g . Let us now consider a typical magnetic field of lpG, 
at 116AU from the star, and vo = 1000 km s~ . We derive 
L = 3 x 10" 9 AU and R = 0.028 AU. Hence L <K R (i.e., 
the Lorentz force dominates the radiation pressure), and p R. 
This yields 2pco g =t 2vo = 2000 km s -1 . The ions turn out to 
have velocities relative to the star randomly distributed over a 
range of 2000 km s~'. Even if the ions do not drift away sig- 
nificantly, their residual velocity range remain far above the 
spectral resolution of the instruments, and they should not be 
observed as a single line at rest with respect to the star. In con- 
sequence, this purely magnetic deceleration model cannot ac- 
count for the observations. 

4.3. Combining gas drag and magnetic field 

We investigate the possible combination of the two preced- 
ing models, i.e., combining a magnetic field and gas drag. The 
equation of motion is now 

dv 

m—=P + qvxB+f , (20) 
df 

where / is the gas drag force. In the general case where / is 
given by Fig. [5] this equation cannot be solved analytically, but 
at low velocity in the induced dipole regime where / = —kv, the 
differential system remains linear and may be solved exactly 
after a somewhat lengthy but straightforward algebra. Keeping 
the definitions and notations of the preceding sections, the so- 
lution is 

v x = cos tb [(Rco g - Lu) r ) cos(wgf) + Lu) g sin(o> g o] e~ Wr ' 

+ Lto r cos tb 

v y = cos ch [Lujg cos(Lo g t) + (Lco r - Rto g ) sin(w f f)J e~ w,f 



Lo) g cos < 



v 7 = sin i 



Roj r ojg — L{co 2 + or) 



+ L sin ( 



Here io g and R are defined as above, co r - lc/m is a fre- 
quency characterising the drag force; L is now defined as 
P = mL(u) 2 + oj 2 ). The motion is still combination of a circu- 
lar motion and a linear drift, but the circular motion is damped 
exponentially at frequency o> r , so that the velocity assumes an 
asymptotic value like in the gas drag case without a magnetic 
field. 

Taking now the numbers given above, and considering a 
typical expected neutral density of 10 4 cm -3 (a column density 



of ~ 10 17 ctrT 2 over ~ 1 AU), we derive to g 
and cog/tor 78. Hence the exponential damping is a slower 
process than the gyromagnetic motion. As co r <K to g , the nu- 
merical value of L is virtually unchanged with respect to the 
preceding definition, and we still have L <k R. After damping, 
the ions have a terminal velocity with respect to the gas given 
by 



L(o4 + co 2 ) 



to 2 + CD 2 sin 2 ( 

\ CO 2 + CO 2 



p 



to 2 + to 2 sin 2 tb 



P . 
* - S1 n< 



(22) 



P/k is the equilibrium velocity v eq given by Eq. ( TTol l in the gas 
drag case, in the case of the induced dipole regime. Hence the 
terminal velocity v t appears in any case less than v eq . It can even 
reach zero if the field is perpendicular to the initial motion. 

Finally, the net result of the interaction is a deceleration 
of the ions that is similar to the non-magnetic case. The ter- 
minal velocity is comparable (or even less), and it is reached 
within the same characteristic time 1 /to r = m/k. The only dif- 
ference concerns the path of the ions. When no magnetic field 
is present, the motion of the ions is linear and they are stopped 
after having encountered a column density of 10 17 cirT 2 , i.e., 
0.7 AU with n = 10 4 ctrT 3 . With a magnetic field, the path 
is no longer linear. When tb — 0, the radial extent of the 
spiral motion of the ions before being stopped is 2p' at 
t — 0, i.e., =i 2R. With the values quoted above, this is about 
0.055 AU; with n = 10 4 crrT 3 , the corresponding column den- 
sity is ~ 8.3 x 10 15 crrT 2 . Hence the ions appear to be stopped 
over a much shorter distance. This is only due to the fact that 
the motion is not linear. As the ions spiral into the quoted dis- 
tance, they encounter the required 10 17 cirT 2 to stop them, but 
not in a linear fashion, rather inside a box of smaller dimen- 
sions. 

The magnetic field appears then as an additional source of 
deceleration for the ions, but the gas drag remains unavoidable. 
Our conclusion is then that a magnetic field may be invoked as 
a refinement to the model, but that it is not necessary, the basic 
process of deceleration remaining the gas drag. 

The main issue concerning the magnetic field should be its 
sudden presence around 100 AU. Why should a magnetic field 
of 1 pG appear there while not present closer to the star ? The 
only possibility is to invoke a kind of heliopause. /3 Pictoris is 
the only normal A type star for which a chromospheric activ- 
ity was detected bv lBouret et alj {2002). These authors derive a 
mass loss rate M = 2.5 X 10 -14 M yr~', with a terminal veloc- 
ity of 200kms~'. Roughly speaking, the suspected heliopause 
should be expected where the magnetic pressure B 2 /2po of the 
surrounding galactic fiel d equals the kinetic pressure of the 
wind. With the values of iBouret etal] d2002h and B = 1 pG, 
it occurs at 529 AU; alternatively, if we want this to occur at 
116AU, a galactic field of 4.5 pG is required. Inside this cav- 
ity, the field would be radial and have no effect on the motion of 
the ions. These are likely values, so that the possibility cannot 
be excluded. 
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5. Conclusion 

The presence of metallic ions at fairly high latitude ov er the 
mid-plane of the f5 Pic circumstellar disk, as observed by B04, 
can be very well explained as a consequence of the FEB pro- 
cess. Whenever the evaporating bodies enter the star grazing 
regime, they are subject to inclination oscillations up to ~ 30 - 
40°. The Can ions released by the FEB during this phase start 
a free, almost radial expansion pushed by a strong radiation 
pressure, keeping track of their initial orbital inclination. Iron 
is also concerned by this process, and we expect Fe i ions to 
be pre sent a t high latitude together with Ca n. The Fe i in the 
data o f lBOl is not null where Ca n is detected. We nevertheless 
explain in our emission line analysis that, unless the electronic 
density is high, the Fe i emission is expected to be weaker than 
the Ca n one, thanks to a predominant ionization of Fe i into 
Fe n. 

This process does not concern Na i ions because, once pro- 
duced by the FEBs, they are quickly photoionized into Nan 
and subsequently no longer experience any noticeable radia- 
tion pressure. Hence we explain the absence of Nai emission 
at high latitude. 

Blown away by strong radiative pressure from the star, the 
Can ions reach the distance of ~ 100 AU in about 1 yr with fi- 
nal velocities of ~ lOOOkms 1 . They need thus to be slowed 
down in order to gather at the star velocity and to form an ob- 
servable line. This can be achieved if the ions encounter at that 
distance a neutral gaseous medium, in agreement with the con- 
clusions of lFQq . A rough estimate of the incoming ion flux due 
to the FEB activity shows that it can account for the necessary 
heating source to render the lines observable. 

In addition to the induced dipole drag force considered by 
lF06l we estimated additional braking effects arising for rapid 
collisional velocities. If we consider the effect of repulsive core 
and inelastic interactions discussed in Sec 4.1.2 and 4.1.3, a 
column density of 10 17 cirT 2 of H i is sufficient to stop the ions 
over a distance of a few AU. The inelastic model we introduced 
is probably very crude, but it is still likely to underestimate 
the effective drag force. Therefore irrespective of the detailed 
description of the interaction processes, the required column 
density r emains below the upper detection l imit of 10 ls cm" 2 
given by Lecavelier des Etangs et al. ( 2001 ). Following lF06l 
it should even be less if we took into account the collective 
plasma behaviour due to partial ionization of the neutral gas 
into account. Conversely, due to the high latitude over the dust 
disk , we do not expect collisions with dust grains (invoked by 
iFOd as a possible braking mechanism) to play a significant role 
in the decelerating process of the incoming ions. 

We also investigate the possible role of a magnetic field 
in stopping the ions. While the sole action of a magnetic field 
is unable to sufficiently slow down the ions, magnetic inter- 
actions provide an additional braking process to the basic gas 
drag model invoked. Combining gas drag and magnetic inter- 
actions can thus be a very efficient way to decelerate the ions, 
still reducing the requirements on the neutral gas density by an 
order of magnitude. This nevertheless constitutes a refinement 
of the model, as gas drag in itself is sufficient to account for 
current observational constraints. 



The key parameter in this model is the distance (~ 100 AU) 
at which the ions are stopped. In the gas drag model, we need to 
assume that no neutral medium is present at 30° inclination up 
to that distance, so that the ions can freely expand radially, and 
that they suddenly encounter some medium there. This would 
mean that the disk begins to significantly flare at that distance. 
As explained above, this distance corresponds also to the ex- 
pected outer edge of the planetes i mal d isk that produces the 
dust, according to lAugereau et al.1 (1200 11) . These two facts are 
probably related. 

Our conclusion is thus that the proposed scenario is plau- 
sible. Another important issue in this study is the number of 
Can or Fe i ions necessary to account for the observations of 
B04. It cannot be determined easily even if we may estimate 
the incoming flux, as it depends highly on the time the ions 
stay within the neutral medium before diffusing away, and sub- 
sequently on the small asymptotic drift velocity they reach. If a 
magnetic field plays a role, this velocity is expected to be sig- 
nificantly lower than without a field; hence the ions should drift 
more slowly across the neutral medium. At a given epoch, for 
the same incoming Ca n flux more ions are therefore expected 
to be trapped in the neutral gas if magnetic forces are active 
than in the opposite case. This is why deriving an incoming 
Ca ii flux in order to compare to the expected number of FEBs 
is very imprecise. This could be the purpose of future inves- 
tigations. Fortunately the uncertainties in the proposed decel- 
eration models are irrelevant here, because once the ions have 
been decelerated, the analytical induced dipole model should 
be valid. 

Our estimate of the incoming ions flux due to FEB activity 
(Sect. 2) is very rough, mainly because the FEB activity itself 
is hard to constrain. Moreover, we expect this activity to be 
time-variable, as c hanges have been observed between various 
observing epochs (Tobi n et"aT]l2004 ). As the emission lines are 
supposed to depend on this flux (via the heating source), we 
expect the strength of the emission lines to present temporal 
variations (at least those at high latitude). It would thus be of 
interest to initiate a follow-up of these lines to check for tem- 
poral changes. 

The FEB scenario is reinforced by the present analysis. The 
off-plane presence of some metallic species, and the absence 
of some others, appear as a natural consequence of the FEB 
scenario and of the mean-motion resonance model with a giant 
planet. This strengthens our view of the /3Pic system as a young 
planetary system. 

Appendix A: Mean-motion resonance theory 

We describe here the non-planar restricted three-body problem, 
with a mass-less test particle orbiting a star, and perturbed by a 
planet orbiting the star. The position vectors of the particle and 
the planet relative to the star are noted r and r' respectively. 
As usual, we call fi the ratio of the mass of the planet to the 
total mass of the system m. We assume that fi <k 1. In this 
framework, the Hamiltonia n of the problem is (see, for instance 
Morbidelli & Moonsll 1993b 



•Ho 
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2a 
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r ■ r 



(A.l) 



Beust and Valiron: High latitude gas in the /? Pictoris system 



15 



where a is the osculating semi-major axis of the orbit of the 
particle. Here the Hamiltonian has been normalised by the con- 
stant factor Qm, where Q is the constant of gravitation. The 
usual way is then to start from the classical canonically conju- 
gate Delaunay variables : 



M , L = 7(1 -fi)a 

cj , G = LVl -e 2 
£2 , H — G cos ; 



(A.2) 



Here M is the mean anomaly along the orbit of the particle, a> 
is its argument of periastron, £2 its longitude of ascending node, 
e its eccentricity and i its inclination with respect to the orbital 
plane of the planet. 

We assume here that the particle is locked in a (p + q) : p 
mean-motion resonance with the planet. It is then of interest 
to introduce new canonically conjugate angle-action variables 
that take into account the resonance : 



cr — A A - m , S — L — G 

q q 

a- = illx -2-A.-Q , S Z = G-H 
q q 



p + q P i 
-v = A A 

q q 



N : 



p + q 



L-H 



The modified Hamiltonian is 



(A.4) 



Here A = M + w + Qis the mean longitude of the particle 
along its orbit; A' is the same for the planet; m — co + Q. 
is th e longitude of peria s tron. T hese variables are i ntroduced 
by iMorbidelli & Moonsl (1 1993b and iMoonsI dl994l) . The an- 
gle cr is often called the "critical angle of the resonance" 
(Moons & M orbidellilll995 ). Resonant orbits are characterised 
by a libration of cr around a stable position. 

The secular dynamics is investigated by performing a 
time-averaging of H over the only remaining fast variable, 
A'. If the orbit of the planet is circular, then the averaged 
Hamiltonian turns out to be inde pendent of v, showing that 
N is a secular constant of motion ( Morbidelli & Moonsl 1993 



Moons & Morbidelli|[l995l) . This can be checked with explicit 
expressions, but this arises from the d'Alembert rules : cr and 
cr z are independent of any axis rotation within the orbital plane 
of the planet, while this is not the case for v. If the planet's 
orbit is circular, the whole Hamiltonian is expected to be in- 
variant for any rotation in that plane, and should consequently 
not depend on v. In that case, the Hamiltonian has two degrees 
of freedom. If we restrict our study to planar motion, then the 
variables cr z and S z disappear and the averaged Hamiltonian 
is integrable. The secular motion is characterised, together 
with the librations of cr, by coupled oscillations in the (a, e) 
plane around the equilibrium value, along a curve N = est. 
This dynamics is de s cribed for many specific resonanc es by 
Morbidelli & Moonsl dl993l) : lMoons & Morbidellil J 1995b . 

If the orbit of the planet is not circular, the action N is no 
longer constant. It is thus able to evolve, but on a much longer 
time-scale than the main librations of cr. Therefore, on a short 
time-scale the oscillations in (a, e) space are preserved, but on a 



longer time-scale, the value of N is subject to change s that may 
drive t he eccentricity to high values. As quoted by I Yoshikawal 
( 1989b . these cha nges are particularly importa nt for resonances 
4:1, 3:1 and 5:2. iBeust & Morbidellil d 1996b showed that the 
4: 1 resonance is a potential so urce of FEBs via this mechanism, 
and Thebault & Beustl ( 2001 ) show that the 3:1 resonance may 
also contribute to the FEB phenomenon. 

If we return to the spatial problem and give an ini- 
tially moderate inclination to the particle, these dynam- 
ics are preserved. In fac t all the simulations present ed in 
Beust & Morbidellil d2000b and lThebault & Beusd d200lb were 
three-dimensional, and the behaviour reported was in perfect 
agreement with the planar description. However, the planar 
model does not describe the inclination oscillations observed 
whenever the eccentricity reaches high values. In order to ex- 
plain them, we must consider the spatial problem as a whole. 



To further study the problem, Morbid elli & Moonsl (1993) 
introduce the following canonically conjugate action-angle 
variables : 



(A.3) ^ 



2tt 1 
To- 2jt 

lf/y = V~ Pvilf/o-, Jo-, J z , Jy) , Jy - N 



S dcr 



(A.5) 



where is the libration period of cr, t is the time, and is 
computed over one libration cycle of cr. p v and p z are functions 
that are periodic with zero average in i/v. 

We are interested in the secular dynamics inside the res- 
onance; hence we perform a second averaging of 'H over if/^. 
The function p v and p z disappear then and J a is a new secu- 
lar constant of motion. Jo- is the normalised area enclosed by 
the libration trajectory in (S, cr) space. It is close to the ampli- 
tude of the libration in cr. Thus in the non-circular problem, the 
value of may change, but the amplitude of the libration is 
roughly pre served. Curves of Jo- = est f or various resonances 
are giv en in Morbidelli & Moons ( 1993 ): Moons & Morbidellil 
dl995b . 

The inclination oscillations are related to the coupled evo- 
lution of tff z and J z . This is not easy to describe in the general 
case as the Hamiltonian fi, even averaged over i/v, is still not 
integrable, as depending on several angles, and because the re- 
lationship between if/ z , v and Jo- and the usual orbital elements 
is not straightforward. A convenient way to investigate the dy- 
namics is to restrict the study to the case J a - 0, i.e., orbits with 
negligible libration amplitude. In this case, the semi-major axis 
a assumes a fixed value (the pericentric equilibrium) close to 
the unperturbed value of the resonance; cr also is fixed to an 
equilibrium value cr . The value of cr depends on the reso- 
nance under consideration. It is then easy to show that under 
these conditions if/, — <r z — cr + co. This method of consid- 



ering zero amplitude libration was used in lBeust & M orbidelli 



( 1996) to draw Hamiltonian maps in the planar problem. 

If we consider the circular (but non-planar) problem, the 
secular Hamiltonian <H, once averaged over if/ IT , has only one 
degree on freedom left. It is a function of if/ z and J z , or alterna- 
tively of i and co, the constant value of N acting as a parameter. 
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Appendix B: The drag force theory 

Given the interaction potential U(r) (the potential energy being 
fj.U{r)), the encounter between a Ca n ion and an H i atom may 
be studied using energy and momentum conservation : 



\ (r 2 + r0 2 ) + U(r) = \v 2 



r z = bv 



(B.l) 
(B.2) 



where the relative motion is described by its polar coordinates 
(r, 0), v is the incoming velocity, b is the impact parameter, and 
r and 6 are as usual shorthands for dr/df and d9/dt. Eliminating 
leads to the radial equation 



r 2 = -2[/(r) + v 2 |l- -| - <,(, ! 



(B.3) 



The classical procedure is then to search for a minimum ap- 
proach distance r m . This is done demanding r 2 = g(r m ) — in 
the preceding equation. There may be no root for r m . Let us for 
example consider the dipole induced potential U{r) = V(r)/fi 
where V(r) is taken from Eq. (0). The resulting equation for 
r m is quadratic, and it is a matter of straightforward algebra 
to show that a root exists only if b > b^, where bo is given 
by Eq. ©. In the opposite case, the distance is expected to 
decrease to zero. However, with a more realistic potential pre- 
senting a repulsive core at small distance, there will always be 
a root for r m : at infinity, Eq. (IB. 3b gives g(r) = v 2 > while at 
small distance, as U(r) > 0, obviously Eq. ( IB. 3b gives g(r) < 0. 

r m corresponds to a polar angle m , and it is easy to see that 
the defection angles is related to 8 m by 



X = 20 m + n 



(B.4) 



m itself can be deduced from the fact that at infinity before the 
encounter, we have = -jr. 

r +o ° de , r + °° e , r +M bv dr m „ 

m +n= \ — dr= 7-dr= ___,(B.5) 

Jr m dr J, \r\ J, r 2 J^T) 



r) 

where g{r) is defined in Eq. dB.3l l. Setting y = r m /r, this ex- 
pression may be rewritten as 



m +Ji 



h - r 



dy 



,(B.6) 



T + 7 b 2 U(r m ) - U(r m /y)] 
from which we derive straightforwardly 

n, ^ 2b f' 1 d >' 

rm Jo ^/l - y 2 + e (b,v,y) 

where 

<*.v,y) = — 

v L 1 - y L 



(B.7) 



(B.8) 



This last expression is of practical interest to numerically com- 
pute x f° r an y potential: e(b, v, y) has a finite limit towards 
y — * 1, so that the presence of the kernel 1/ yl - y 2 allows a 
rapid computation using a Gauss-Chebychev quadrature tech- 
nique. 



Once x is known, the impulsion change to the ion during 
the encounter reads 



dp = fiv 



cosx — 1 
sin^cos <f 
sinx sin (p 



(B.9) 



once written in a Cartesian (x, y, z) referential frame with the 
x-axis parallel to the initial motion of the ion, and where <p is 
an azimuthal angle characterising the plane of the encounter. 
Of course over many encounters is a random angle, so that 
the y and z components of dp vanish. On average then, we have 



dp = n(cosx - 1) v 



-2,sin 2 (f) 



(B.10) 



Then follows a classical cross-section calculation. The number 
of encounters occurring within the time-span 6t at (b, (f>) within 
db and d(f> is nvdtbdbdO. Multiplying then dp by this number 
of encounters and integrating over and b, we derive the drag 
force as: 



/ 



hi! 



5p(b,v)nv6tbdb d0 



— —4nn/jv 



Jo 



b sin — - — db 



(B.ll) 



If the collision is not elastic, the modulus of the relative veloc- 
ity v' after the encounter differs from the initial one v. We have 
v' = ev. The impulsion change to the ion during the encounter 
reads now 



8p = n 



vcos^ - ev 
vsin^cos^ 
v sin^f sin (f> 



(B.12) 



As in the elastic case, there is no average change perpendicular 
to the motion. In the direction of the motion, we derive 

ty = (-2// e sin 2 (!) + ( e -l)//)v . (B.13) 

Integrating as above over b and we derive the drag force as 



/ = -Inn^iv 



2c sin 2 I 1 + i - c \bdb 



.(B.14) 
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